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Abstract. We analyze several aspects of the phenomenon of stochastic resonance 
in reaction-diffusion systems, exploiting the nonequilibrium potential's frame- 
work. The generalization of this formalism (sketched in the appendix) to extended 
systems is first carried out in the context of a simplified scalar model, for which sta- 
tionary patterns can be found analytically. We first show how system-size stochas- 
tic resonance arises naturally in this framework, and then how the phenomenon 
of array- enhanced stochastic resonance can be further enhanced by letting the 
diffusion coefficient depend on the field. A yet less trivial generalization is ex- 
emplified by a stylized version of the FitzHugh-Nagumo system, a paradigm of 
the activator-inhibitor class. After discussing for this system the second aspect 
enumerated above, we derive from it — through an adiabatic-like elimination of 
the inhibitor field — an effective scalar model that includes a nonlocal contribu- 
tion. Studying the role played by the range of the nonlocal kernel and its effect 
on stochastic resonance, we find an optimal range that maximizes the system's 
response. 



1 Introduction 



Stochastic resonance (SR) is nowadays a paradigm of the constructive effects of fluctuations on 
nonlinear systems [1)2] . Sketchily, the phenomenon occurs whenever the Kramers' rate for the 
transition between attractors matches the typical frequency of a signal which is incapable by 
itself to trigger that transition (i.e. it is subthreshold). Whereas several measures of SR can be 
defined [the signal-to-noise ratio (SNR) and the spectral amplification factor (SAF) being the 
main ones], theoretical analysis is usually carried on in terms of the two-state approximation 
[T]. Since its discovery a quarter of century ago — and besides exploring related phenomena 
like e.g., coherence resonance |3] — interest has gradually shifted towards increasingly complex 
systems, networks and nonlinear media being the main directions. Instances of this trend are 
the experiments carried out to explore the role of SR in sensory and other biological functions 
and experiments in chemical systems [5]. 

Our concern throughout this review will be with nonlinear media that can be described 
as reaction-diffusion (RD) systems, namely those that can be thought of as a collection of 
diffusively coupled nonlinear units. The possibility of enhancing the system's response through 
the coupling of those units [6171819110111112] has been among the issues explored during the 
last decade, together with the "naturalness" problem (how does nature manage to make the 
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system's response less dependent on a fine tuning of the noise intensity) or that of searching 
for different ways to control the phenomenon |13ll4j . 

In dissipative dynamical systems, the very notion of Lyapunov's function is as useful as 
that of attractor itself. Even when often it cannot be explicitly computed (because integrability 
conditions are not readily met), it allows for picturesque reasoning in terms of "energy land- 
scapes" or "attraction basins" . When those dynamical systems are submitted to forces that can 
be modeled (a la Langevin) as stochastic, a new meaning — statistical in nature, akin to the 
notion of "free energy" — is added to the picture (in fact, it is worth mentioning that the first 
function known to have the Lyapunov property was Boltzmann's H- function) . Moreover, even 
for vanishing noise intensity, the very existence of stochastic terms (a "transport matrix" ) can 
render the system well conditioned regarding integrability conditions. That was the rationale 
behind the definition of nonequilibrium potential (NEP) [16 17 18 19 , two approaches to which 
are described in the Appendix. Such NEP is a special Lyapunov's function of the associated 
deterministic system, which for nonequilibrium systems plays a role similar to that played by 
a thermodynamic potential in equilibrium thermodynamics |16j . It is closely related to the sta- 
tionary solution of the system's Fokker-Planck equation, and characterizes the global properties 
of the dynamics: attractors, linear and relative stability of these attractors, height of the barri- 
ers separating attraction basins. In addition, it allows to evaluate the transition rates among 
the different attractors [16117118119120] . Regarding the problem of SR in extended systems, it 
was shown that the knowledge of the NEP allows to obtain a rather complete picture of the 
behavior of the output signal-to-noise ratio (SNR). The novelty with nonequilibrium extended 
systems is that even pointlike attractors in the medium's infinite-dimensional phase-space can 
be nontrivial field configurations (real-space patterns). 

In a series of recent papers we have studied the SR phenomenon for the transitions between 
two different patterns |8l9ll0lllll2ll5j . exploiting the concept of nonequilibrium potential. In 
this review we discuss some recent results concerning different aspects of SR in RD systems. In 
Sec. [2] we discuss the phenomenon of system-size stochastic resonance (SSSR), and show how 
can it be analyzed and understood within a NEP framework |21j . In Sec. [3] after reviewing a 
recent study on the enhancement of the SNR found for a scalar system with density-dependent 
diffusivity [12], we discuss its extension [22 to an array of FitzHugh-Nagumo units [53]. In 
Sec. [i] — through an adiabatic-like elimination of the inhibitor field in an activator-inhibitor 
system — an effective scalar system with a nonlocal term is derived, and the role of the local and 
nonlocal interactions on the SR response studied. The main conclusions are finally summarized 
in Sec. |5] 

2 System-size stochastic resonance 

Recent studies on biological models of the Hodgkin-Huxley type 24 25j have shown that ion 
concentrations along cell membranes display intrinsic SR-likc phenomena as the number of ion 
channels is varied. A related result [25] shows that even in the absence of external forcing, the 
regularity of the collective firing of a set of coupled excitable FitzHugh-Nagumo units is optimal 
for a given number of elements. From a physics point of view, the same phenomenon — called 
system-size stochastic resonance (SSSR) — has also been found in an Ising model as well as in 
a set of globally coupled units described by a 4> 4 theory [27] . It has been even shown to arise in 
opinion formation models |28j . 

Since the SSSR phenomenon is peculiar to extended systems, there is an obvious interest 
in describing it within a NEP framework, that offers a very general framework for the study 
of the dependence of SR and related phenomena on any of the system's parameters. Here we 
discuss in some detail a one-component ("scalar") RD system [2T], and briefly refer to other 
cases analyzed in [27] and [29] . 

2.1 Review of the scalar model 

For one-variable dynamical systems, the Lyapunov's function can always be found by quadra- 
ture. This property can be readily translated to scalar RD systems: the Lyapunov's functional 
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J-[4>] fulfills the "potential" condition d t <j)(y,t) = — SJ r [(f)]/S(j)(y, t), where 5/Scj)(y,t) indicates a 
functional derivative. This is also the NEP for a scalar transport matrix (i.e. a multiple of the 
unit matrix in the medium's infinite-dimensional phase space). 

The specific model we shall focus on here has a piecewise linear reaction term, that mimics 
general bistable RD models [23], e.g. those with a cubic-like reaction term. In the following, 
we shall exploit some of the results on the influence of general (partially reflective or albedo) 
boundary conditions found in [30., as well as previous studies of the NEP [T7] and of SR 
[819111115] . The particular dimensionless form of the deterministic model we start with is [301819] 



dt 



D 



(Pi 
dy 2 



(1) 



where 0(z) is the Heaviside step function, 4> c is the value at which the piecewise-linear "reaction 
term" has t he ju mp, and D is a phenomenol ogic al "diffusion coefficient", not necessarily related 
to 7 in Sec. 



2.2 



[in Graham's approach, Eq. (15), D/(Ay) 2 would be the matrix elements Q vv±l 
in a discretization of the Laplacian] . All the effects of the parameters that keep the system away 
from equilibrium (such as the electric current in electrothermal devices like the ballast resistor 
23 30], or some external reactant concentration in chemical models) are now included in <p c . 
For the system to display a bistable behavior, it must be < cp c < <fih- 

We consider here the class of static structures <fi(y) studied in [3U]. They are even solutions 
to the stationary (d t <fi — 0) version of Eq. (HI in the bounded domain [— j/Lj^lL with equal 



albedo boundary conditions (b.c.) at both ends 

d<Ky,t) 



dy 



= Tk4>(±y L ,t). 

v=±vl 

k > is called the albedo parameter: the limit k — > yicds Neumann's b.c. dy<fi(y,t)\ y —± yii = 
and the k — > oo one, Dirichlet's b.c. 0(±y/,,t) = 0. 



The explicit form of these static patterns is 

sinh^y (k,*$) //>(*, ;J%) 
,.>[,,) = < >,, { I - cosh(y)p (k, ^=g£) lp (k, $l 



sinh(y c )p> (fc,M)/p(fc,-Mfc) 



-yL <y < -y c 
-y c <y < y c , 
y c <y < vl, 



(2) 



where p(k,() = sinh(^) + k cosh(£) and p'{k,Q — dp/d(. The coordinate values yf at which 

4>{Vc) = <i>c are 
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<\>h 



p k, 
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Each real solution yf < y^ to Eq. ([3| represents a structure with a central "activated" zone 
(4> > 4> c ) and two lateral "resting" regions (0 < <fr c ). Figure 5 in [30 displays the relation y c /yL 
vs k, for several values of 4> c /4>h. 

Typical shapes of the arising patterns are shown in Fig.[T] Through a linear stability analysis 
it has been shown [30] that the structure with the smallest "excited" region [that is with 
y c = y+, denoted by 4> u {y)] is unstable, whereas the other one [with y c — y~ , denoted by <fii(y)] 
is linearly stable. The trivial homogeneous solution 4>o(y) — exists and is linearly stable for 
any parameter set. These two linearly stable solutions (<j>o and (pi) are the only stable static 
structures under albedo b.c. We will concentrate on the region of values of (j> c /4>h, Vl and k 
where <pi exists. 

For the finite system with albedo b.c, the NEP is a functional of 4> and a function of k, y^ 
and (pc/'Ph- It has the expression [T7] 



T{[<j)},(j) c /(j> h ,k,y L ) 





r<f>(y,t 
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{-ft + <p h 0(0' - cj> c )] + 
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Fig. 1. Inhomogeneous static solutions to Eq. for y^ = D = 4>h = 1.0 and <f) c — 0.193 (highlighted 
by the dashed horizontal line). Dash-dotted upper curve: 0i(y) for k = 1.0 (in this case, y c = y^ > 
Lower curves: </>i(y) (solid line, for which +y~ is highlighted) and (f) u {y) (dotted line, for which — 1/+ is 
highlighted) for k — 7.0. 




Fig. 2. NEP J-([<j)], k, j/i) evaluated at the inhomogeneous stationary solutions cf>i{y) (lower branch) 
and 4> u {y) (upper branch), as a function of: (a) system's size 2/l, with k = 3.0; (b) albedo parameter k, 
with yL = 1.2. The remaining parameters are D = 1.0, 4> c /4>h = 0.193. The NEP for the homogeneous 
stationary solution 4>o(y) coincides with the horizontal axis. 



When the NEP is evaluated at the inhomogeneous static solutions of Eq. ([I]) [Eqs. ^ and 
^] it takes the explicit form |8I17] 



F uA (<j) c /(j) h ,k,y L ) = F{[<i)u,i],<t>cl '<Ph,k,yL) 

2(j) c 



(4) 



1 - 



sinh [yf/VD 



P (k,y L /VD 



while at the trivial solution <j) = it is !F([cf)o], Vi) = = 0. 

Figure [2^l depicts the nonequilibrium potential F([<f>],yL) as a function of the system's size 
yL, keeping the albedo parameter k and the ratio <p c /<Ph fixed. The curves correspond to the NEP 
^^(Vl), whereas T° coincides with the x-axis. Our focus is the "bistable zone" y L ^ 0.72, 
where 4>\{y) exists. The unstable structure 4> u (y) is a saddle point for F[<fi\ (in the medium's 
infinite-dimensional phase space), so its NEP T u (%il) > (upper branch in Fig. |2k). On the 
other hand, both </>i (y) (lower branch) and 4>o (iB-axis) are local minima of the NEP. 

One immediately notices that AT l {yL) = ^F u {yi,)~ F \Vl) is an (almost linearly) increasing 
function of y^ (this has a profound implication for SSSR, as we shall see below) . Equation ([3]) has 
real solutions only for y^ ^ 0.72. This corresponds to a supercritical saddle-node bifurcation, 
at which both inhomogeneous structures pop up. Now, the most important feature in Fig. |2^i is 
that J- (yr,) vanishes at a certain system's size y* L {~ 1.0 for the given values of k and <fi c /4>h)- 
At that point, the stable inhomogeneous structure 4>\{y) and the trivial solution 4>o(y) exchange 
their relative stabilities. 
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For completeness, in Fig. [2]d we plot J-([<f>],k) for the same value of 4> c /4>h as in Fig. [2^. 
For the chosen value of it is always ^ rl (fc) < 0, and correspondingly there is no "stability 
exchange" as a function of k. Also, the initially large AJ rl (k) decreases with k as ^-""(fc) and 
J rl (fc) tend (for k — > oo) to the values corresponding to Dirichlet's b.c. [17. 



2.2 Results for SSSR 

By including an additive spatiotemporal noise source £(y, t) |15l33j . Eq. becomes a stochastic 
partial differential equation for the random field <$>(y, t). The simplest assumptions about t) 
are that it is Gaussian, with zero mean and a correlation function given by t)£(y', f')) = 
2 r y8{t — t')S(y — y'), where 7 denotes the noise strength. 

As discussed in |8|9|11|15] . known results for activation processes in multidimensional sys- 
tems [3JJ allow us to estimate the activation rate using the following Kramers'-like expression 
for the mean first-passage time for the transitions between attractors 



<Tj) = r exp 



where AJ^{y£) — J- u (yL) — ^(vl), i = 0, 1. The prefactor To is usually determined by the 
curvature of !F\<j)\ at its extrema. On one hand, it is typically several orders of magnitude smaller 
than the average time (t) , while on the other it does not change significativcly when varying the 
system's parameters around the "bistable point" yj, where .F([<^o], y*h) = 3~([4>\\, y* L ). Hence we 
may simplify the analysis by assuming here that To is constant, and scale it out of our results. 
The behavior of (t) as a function of k and 4> c /4>h has been shown in |8l9ll7j . 

As done in [5], we now assume that the system is (adiabatically) subject to an external 
harmonic variation of the parameter <f> c : <j> c (t) = </>° + Scj) c cos(wt) |9ll5j . and exploit the "two- 
state approximation" [1] as in [9 11 15]. Such approximation reduces the whole dynamics on 
the bistable potential landscape to one where the transitions occur only between the states 
associated to the bottom of each well, hence the only dynamical contents resides in the transition 
rates. Up to first order in the amplitude 5<fi c (assumed to be small in order that the periodic 
input be sub-threshold) the transition rates Wj adopt the form 



w ,4 



Hi =F cti — - cos(tut) 
7 



(5) 



where (at constant </>/,) fXi w exp[— AJ 7% (<j)®,yL)/'f[ and on « fii(dAJ 7% /d^) c \^o), i — 0,1. The 
quantity inside parentheses can be obtained analytically using Eq. Q. These results allow to 
calculate the autocorrelation function, the power spectrum density and finally the SNR, that 
we indicate by R. The detailed calculation can be found in the appendix of [TT]. Up to the 
relevant order (the second) in the signal amplitude 8(j) c , we obtain 

tt (a Mi + "iMo) 2 TT ^! 

R = ~, = ^> (6) 

4/xo/xa ^0 + Mi ^1 + Mi 

where we have used the form of the on to reduce the expression, and defined <P = \24>h 2/ c (?/l)] 2 - 
Figure [3] (left) is a plot of R as a function of the noise intensity 7 for a fixed system's length 
y^, displaying the typical maximum that has become the fingerprint of the SR phenomenon. 
In Fig. [3] (right), the roles of 7 and yi, are exchanged (R is plotted as a function of yL for fixed 
7). Such a response is the expected one for a system exhibiting SSSR. In both cases, the values 
of k and 4>®/4>h are kept fixed. 

Within the NEP context and in this kind of systems, the phenomenon arises due to the 
breakdown of the NEP's symmetry. This means that (as shown in Fig. SI when varying y^, 
both attractors can exchange their relative stability. For y^ = y* L ~ 1 both stable structures — 
the inhomogeneous one <f>i(y) and the trivial one 4>o — have the same value for the NEP. For 
Vl < y*Li 4>i{y) becomes less stable than 4>q so transitions from <f>i(y) to 4>o are more frequent 
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Fig. 3. Left: SNR vs the noise intensity 7, with system's size i/l ~ 1.1. Right: SNR vs i/l, for 7 = 0.1. 
The remaining parameters are k — 3.0, D = 1.0, and <j>i/4>h = 0.193. 



(the barrier is lower) than in the reverse direction, thus reducing the system's response. When 
Ul ~ 0.72, 4>i(y) and 4> u {y) coalesce and disappear, and the response is strictly zero (within 
the linear response scheme implicit in the two-state approximation). When yi > y* L , 4>i(y) 
becomes more stable than 4>o, making now transitions from <po to <fii(y) more frequent than in 
the reverse direction, and reducing again the system's response. Clearly, the system's response 
has a maximum when both attractors have the same stability (j/l = y* L ), and decays when 
departing from that situation. Hence, for this system and within this framework, SSSR arises 
as a particular case of the more general discussion done in |llj . It should not come as a surprise 
to find an analogy with the mechanism of double stochastic coherence described in |32j , where 
the NEP's symmetry is induced by (an additional, multiplicative) noise. 

By comparing figures^ and [3] it becomes apparent that the value of y^ at which the SNR 
has its maximum differs slightly from y* L (where the crossing between T 1 and T° takes place). 
The origin of this discrepancy is the following: whereas on qualitative grounds we have argued 
that the maximum of the SNR should be related to the potential being symmetric (both wells 
having the same "energy") [TT], the exact condition is that the transition rates between both 
wells be equal. In general, due to small differences between the curvatures at the bottom of each 
well, those rates become equal for values of yi, slightly different from the one at the symmetric 
case. Although by adopting here a constant value of To we have assumed equal curvatures, there 
is still a difference in the values of the on, since the dAT r / d<t> c j^o, i = 0, 1 are slightly different 
(a fact reflected in the dependence of <P on y^). 

Additional light can be shed on the phenomenon when viewed from a different angle. Figure 
[4]is a plot of R as a function of k at fixed values of 7, yj, and cfP c / cf>h. It exhibits a broad resonance 
since for k not too large (indicating a high reflectiveness at the boundary or a reduced exchange 
with the environment) R increases with k, whereas it slightly decreases for larger k values (where 
the system's boundaries become absorbent). An explanation of this behavior in terms of the 
NEP has been given in [29] : as already observed, the NEP's symmetry is broken for this value of 
yL', moreover, whereas the lower branch in Fig.[2]D goes rapidly towards the value corresponding 
to Dirichlet's b.c. the upper branch keeps increasing, thus degrading the SNR. In any case, the 
fact that the resonance is broad indicates the robustness of the system's response with regard 
to fc, a parameter that (together with 7) encodes the coupling with the environment. 

We stress the fact that the NEP framework put forward in this review allows to study 
SSSR between whole patterns. The explanation offered in |27] to the phenomenon resorted to a 
collective variable X w (1/-/V) S J= i x ji an d to the fact that the noise in the effective stochastic 
differential equation for X scaled with size. In it was shown that all the cases discussed in 
[2"T] can be put within the same NEP framework than the above studied scalar model. In fact, 
the aforementioned almost linear increasing dependence of AJ- 1 (yr / ) on yj, can be interpreted 
as a noise scaling with size. There are however situations where the NEP's symmetry is retained 
as the system's size is varied. We may then speak of a genuinely noise-scaling SSSR, in contrast 
to the cases that could be called NEP symmetry breaking SSSR [29]. 
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k 

Fig. 4. SNR vs k for y L = 1.2, 7 = 0.1 and (f>°/4> h = 0.193. 

3 Case of selective coupling 

In this section we analyze SR in two extended systems with density-dependent diffusive- like cou- 
pling: an extension of the scalar RD model considered in Sec. [2] [12J, and an array of FitzHugh- 
Nagumo [23] units. 



3.1 Scalar model 

Here we extend the one-component RD model discussed in Sec. [2] by letting the diffusive pa- 
rameter D in Eq. ([I]) depend on the field 4>{x, t). As a matter of fact, since in the ballast resistor 
|23l30j the thermal conductivity is a function of the energy density, the resulting equation for 
the temperature field includes a temperature-dependent diffusion coefficient in a natural way. 
The form of the governing equation is now 

d t <P(x, t) = d x {D(<f>)d x cf>} + /(<(>) + ^(a:, t), (7) 

with t) and f((f>) as in Sec. [2] 

As it was done for the reaction term, a simple choice (that retains however the qualitative 
features of the system) is to consider the following dependence of the diffusion term on the field 
variable 

D{4,)=D [l + he{(j)- ( j> c )). 

For simplicity, here we choose the same threshold <p c for the reaction term and the diffusion 
coefficient. 

We assume the system to lie in a bounded domain [-L, L), with Dirichlet b.c. at both ends: 
4>(±L,t) — 0. The form of the patterns is analogous to what has been obtained in Sec. [2j the 
only difference being that in the present case d(f>/dx\ Xc is discontinuous and the area of the 
"activated" central zone depends on h. 

As before, the indicated patterns are extrema of the NEP: the unstable pattern <j> u (x) is 
a saddle-point of this functional, separating the attractors 4>o(x) and <fi s (x)- For the case of a 
field-dependent diffusion coefficient D(<p{x,t)) as described by Eq. Q, the NEP reads [T2] 



™-£{-jT 



D{4>')f{ct>')d<t>'+\ 



D ^fx 




Given that dt<j> = —[l/D(<p)](5!F/5(f>), one finds dT jdt = — J (S!F/S(j)) 2 dx < 0, thus warranting 
the Lyapunov's functional property. 

Whereas in Sec. [2] we kept constant and varied yL, we now vary instead <fP c at constant 
L. Similarly as before, both linearly stable states have the same value of the NEP (i.e., they 
are equally stable) at some value <#j_pf the threshold. The way !F[4>] depends on (ffi c resembles 
the dependence on y^ shown in Sec.p] but now AT S is an (almost linearly) decreasing function 
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Fig. 5. Left: SNR vs the noise intensity 7, for Do — 1.0 and h = 0.0 (full line), —0.25 (dashed line) 
and 0.25 (dotted line). Right: the maximum i? ma x of the SNR curve as a function of the selectiveness 
h of the coupling, for Do = 0.9 (dashed line), 1.0 (full line) and 1.1 (dotted line). The arrows a and b 
indicate the response gain due respectively to a homogeneous increase of the coupling and to a selective 
one. The larger gain in the second case is apparent. The inset shows the dependence of _R m ax on Do for 
h — —0.25 (lower line), 0.0 and 0.25 (upper line). The remaining parameters are L = 1.0, 8<f) c = 0.01 
and 12 = 0.01. 

of (jP c and both inhomogeneous structures coalesce and disappear through a subcritical saddle- 
node bifurcation. As in the previous case, we analyze only the neighborhood of 4>* c . We shall 
moreover consider only the neighborhood of h — 0, where the main trends of the effect can be 
captured. 

Figure [5] (left) depicts the dependence of R on the noise intensity 7 for three values of h, 
each curve displaying the typical SR maximum. Figure [5] (right) is a plot of the value -R max 
of these maxima as a function of h. The dramatic increase of i? max (several dB for a small 
positive variation of h) is apparent, and shows the strong effect that the selective coupling (or 
field-dependent diffusivity) has on the system's response. 

It must be noted that the only two approximations made in order to derive our results — 
namely the Kramers-like expression and the two-level approximation used for the evaluation of 
the correlation function |12j — break down for large positive values of h because for increasing 
selectivity the curves of T\<j)\ vs (f)® shift towards the left [T2], which in turn means that the 
barrier separating the attractors at <j>* c tends to zero. This effect is basically the same as the one 
discussed in Refs. [9111] in connection with global diffusivity Dq. It is also worth noting that 
except for the two aforementioned approximations, all the previous results (e.g. the profiles 
of the stationary patterns and the corresponding values of the noncquilibrium potential) are 
analytically exact. 

3.2 FitzHugh-ISIagumo model 

Here we study an array of FitzHugh-Nagumo |23j units, with a density-dependent (diffusive- 
like) coupling. The NEP for this system was found within the excitable regime and for particular 
values of the coupling strength [19]. In the general case, however, the form of the NEP has not 
been found yet. Hence, we have resorted to a study based on numerical simulations, analyzing 
the influence of different parameters on the system's response. Nevertheless, the idea of the 
existence of such a NEP has always underlied this study. The results show that the enhancement 
of the SNR found for the scalar system [TJ] is robust, and that the indicated non-homogeneous 
coupling could clearly contribute to enhance the SR phenomenon in more general situations. 

We consider a simplified version of the FitzHugh-Nagumo model |12|17)23] . which has been 
useful for gaining qualitative insight into the excitable and oscillatory dynamics in neural and 
chemical systems. It consist of two variables: 

— a (fast) activator field u, that in the case of neural systems represents the voltage variable, 
while in chemical systems represents the concentration of an autocatalytic species. 
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— an inhibitor field v, associated (within a neural context) to the concentration of potassium 
ions in the medium, and that in a general chemical reaction inhibits the generation of the 
u species. 

Instead of considering the usual cubic-like nonlinear form, we use a piecewise linear version 



du{x,t) d 

dt dx 

dv(x,t) d 

dt dx 



+ f(u)-v + £(x,t) (8) 
+ /3u-av, (9) 



where f(u) = —u + 0{u — <f> c ), and £(x,t) is a <5-correlated white Gaussian noise, as before. 
7 indicates the noise intensity and <p c is the "discontinuity" point, at which the piecewise 
linearized function f(u) presents a jump. The parameter e = t u /t v -C 1 indicates the timescale 
ratio between the (fast) activator and the inhibitor. Together with a and (3, it is chosen to 
correspond to the excitable regime. We consider Dirichlet b.c. at x = ±L. Although the results 
are qualitatively the same as those that could appear considering the usual FitzHugh Nagumo 
equations, this simplified version allows us to compare directly with the previous analytical 
results for this system |12ll5j . 

As in [12], we assume that the diffusion coefficient D u (u) is not constant, but depends on 
the field u according to D u (u) = [1 + hO(u — <f> c )]. This form implies that the value of 
D u (u) depends "selectively" on whether the field u fulfills u > <fi c or u < <p c . D® is the value of 
the diffusion constant without such "selective" term, and h indicates the size of the difference 
between the diffusion constants in both regions [if h = then D u (u) — constant]. D v (v) is 
the diffusion for the inhibitor v, that we assume to be homogeneously constant. 

This system is known to exhibit two stable stationary patterns. One of them is u(x) = 0, 
v(x) — 0, while the other is one with nonzero values. Further, we consider, as before, that 
an external, periodic, signal enters into the system through the value of the threshold <j> c , 
4> c {t) = 4>° c + 54> cos(u>t), where to is the signal frequency, and S(f> its intensity. 

All the results were obtained through numerical simulations of the system. The continuous 
version of the system indicated by Eqs. (J8j> , ^ was transformed into a second-order spatially 
discrete one 

Ui = D Ut Aui + (D Ui+1 - D Ui _ 1 )(u i+1 + Ui-t) + f(Ui) -Vi + gi(t), 
Vi = D v Avi + (3iii- avi, 

with A(f>i = (4>i—i + 4>i+i — 2cf>i). The extensive numerical simulations performed for a set of 
equations were done exploiting Hcun's algorithm 33J . 

In this spatially-extended system there are different ways of measuring the overall system's 
response to the external signal. In particular, we show the evaluated output SNR in two different 
ways (the units being given in dB): 

— SNR for the middle element of the chain evaluated over the dynamical evolution of ujv/2i 
that we call SNR 2 (however, having Dirichlet b.c, the local response depends on the distance 
to the boundaries). 

— In order to measure the overall response of the system to the external signal, we computed 
the SNR as follows: We digitized the system's dynamics to a dichotomic process s(t): At 
time t the system has an associated value of s(t) = 1 (0) if the Hilbert distance to pattern 
1 (0) is lower than to the other pattern. Stated in mathematical terms, we computed the 
distance T>i\-, ■] defined by 

s. 1/2 



V 2 [f,g] = y 



dx [/(x) - g{x)f 



in C 2 ([— L,L]), the Hilbert space of the real- valued functions in that interval. At time t, 
digitized process is computed by means of 

_ f 1 if V 2 [P?(x),u(x,t)] < V 2 [Ptf(x),u(x,t)] 
sw ~ 1 if V 2 [P?(x),u(x,t)} > V 2 [P "(x),it(x,f)] ' 
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Model 


Numerical 


a /3 e cfp c 8<t> lo D v 
0.3 0.4 0.03 0.52 0.4 5tt/8 1.0 


At N 
IO -3 51 



Table 1. Fixed parameters for the FitzHugh-Nagumo model with D u (u) 




Id' 2 10"' 10° io 1 1<T 10"' io° io' 

7 7 



Fig. 6. SNR 2 and SNR P vs the noise intensity 7, for D° u = 0.3 and h = -2.0 (+), -1.0 (A), 0.0 (o), 
1.0 (□), and 2.0 (O)- Both measures reveal a systematic enhancement of the SNR as h increases. 




Fig. 7. SNR 2 and SNR P vs the selectiveness of coupling h, for D° = 0.3 and 7 = 0.032 (O), 0.32 (□), 
0.6 (o), 1.2 (A), and 3.2 (x). 

We call this "global- like" measure SNR p . 

The parameters kept fixed have been summarized in Table [l] The simulation was repeated 250 
times for each parameter set, and the SNR was computed by recourse of the average power 
spectral density. 

Figure [6] depicts the results for the different SNR measures we have previously defined, as 
functions of the noise intensity 7. For both measures it is apparent that there is an enhancement 
of the response for h > 0, when compared with the h = case, while for h < the response is 
smaller. 

In Fig. [7] we show again the response's measures, but now as functions of h. We have plotted 
the maximum of each SNR curve for D° = 0.3, and 7 = 0.01, 0.1, and 0.3. It is clear that there 
exists an optimal value of 7 for which the response is largest. The rapid fall in the response for 
h < is also apparent. 

In Fig. [8] we show the dependance of SNR on h, for different values of the diffusion which 
depends on the activator density D^. It is apparent that the response becomes larger when 
the value of D° is larger. However, as was discussed in [8 15 , it is clear that for still larger 
values of Z?°, the symmetry of the underlying potential (that is the relative stability between 
the attractors) is broken and the response finally falls down. 

The previous figures clearly show that the response to the external signal grows with the 
"selectiveness" of the coupling, showing the robustness of the phenomenon presented in [12115] . 
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Fig. 8. SNR 2 and SNR P vs h for 7 = 0.032 
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and D° = 0.0 (Q), 0.1 (□), and 0.3 (o). 



4 Nonlocal Interaction 



Let us consider again a system like the one described by Eqs. but now assume that 

D u and D v are constant. In Ref. [18] it was assumed that the inhibitor-like field has a diffusive 
transport behavior, and is fast enough that can be adiabatically eliminated, thus yielding an 
effective scalar RD equation with a nonlocal term, characterized by a diffusive kernel G(x,x'). 
After briefly reviewing the derivation of the NEP for that situation, we shall assume in this 
section that the transport mechanism of the adiabatically eliminated inhibitor-like field is of 
nondiffusive character, thus yielding a kernel H{x 1 x') that is more localized in space, and with 
a controllable interaction range. 

Following Ref. [T5] , let the system be defined by 

dtu(x, t) = D u d 2 u(x, t) — u(x, t) + 0[u(x, t) — a] — v(x, t), 
e _1 dtv(x, t) = D v d x v(x, t) + (3u(x, t) — av(x, t), (10) 

where e was defined after Eqs. ([8]), (|9J, and let it be confined to the domain [—L,L], with 
Dirichlet b.c. at both ends: u(±L,t) = v(±L,t) — 0. Contrarily to the standard hypothesis, 
we now assume that the inhibitor is much faster than the activator (i.e. t v <ti t u ). In the limit 



e — > oo, we can rewrite Eq. ( 10 1 as 



dtu[x, t) — D u d^u(x, t) — u(x, t) + 0[u(x, t) — a] — v(x, t), 
= D v d^.v(x, t) + 0u(x, t) — av(x, t). 

In the last pair of equations we can eliminate the inhibitor (which is now slaved to the activator) 
by solving the second equation using the Green's function method 

[-D v d 2 x + a] G(x, x') = S(x - x'), 

v(x) = (3 I dx'G(x,x')u(x'), 



where the Green's function G(x, x') is given by 

„, ,.1 J [sinh k(L — x')/ sinh 2kL] sinh k(L + x) x<x', 
^' X DJ: { [sinh k(L + x') / sinh 2kL] sinh k(L - x) x > x' , 

with k = (a/Dy) 1 / 2 . This slaving procedure reduces our system to a nonlocal equation for the 
activator only, that has the form 

du(x,t) d 2 u(x,t) f L , 

~ D u — \-J\u)—p I &{x,x )u(x jdx . (11) 



at u dx 2 



L 
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From this equation, and taking into account the symmetry of the Green's function G{x, x'), we 
can obtain the Lyapunov functional for this system, which has the form 



dx 



-L 



Du 
2 



du 
dx 



f(w) dw - 







dx' G{x, x) u(x') u{x) 



-L 



(12) 



This spatial nonlocal term in the NEP takes into account the repulsion between activated 
zones. When two activated zones come near each other, the exponential tails of the inhibitor 
concentration overlap, increasing its concentration between both activated zones and creating 
an effective repulsion between them. Hence the Green's function plays the role of an exponential 
screening between the activated zones. In Ref. [9] the knowledge of such NEP was exploited to 
study SR on the system indicated by Eq. p"T| ). 

The starting point of our present analysis will be the effective, nonlocal and stochastic RD 
equation for the real (activator-like) field 4>(x,t), analogous to Eq. (Ill, defined in the one 
dimensional domain x G [— L, L] by 



50 

at 



= D 



dx 2 



+ /(0) - p / H(x, x') 0(x') dx' + £(x, t), 



— L 



where the diffusivity D is constant, and we assume a cubic nonlinear term /(</>) = (4>—b) (2—0). 
Here £(x, t) is an additive Gaussian white noise, as in the previous cases. As before, the system 
is subject to Dirichlct b.c. <fi(±L,t) = 0. 

Similarly to Eq. ( |11[ ), this system could be written in a variational form, with the func- 
tional T[4>] given by Eq. (12 1. As anticipated, here we consider a nondiffusive kernel, with a 
controllable interaction range. In order to keep our analysis simple we propose the following 
form 

1/2; \x-x'\ < I 
0; | a; - x'\ > I, 



H(x,x') 



(13) 



which allows the analysis by just varying the interaction range 21. 

The new effective RD equation contains local and nonlocal couplings (corresponding to the 
diffusive and the nonlocal contribution, respectively). The last one contains the nonlocal kernel 
given by Eq. ( |T3| , with a variable range 21, that corresponds to the interaction of the field at 
points x' € [x — l,x + I]. However, such points will contribute if and only if they are inside 
the domain [— L, L]. We are now in position to study the role played by the nonlocal kernel 
(particularly by its range 21) on the SR phenomenon. 

As before, the SR between stationary solutions was investigated in terms of the two-state 
approach (all the details about the procedure and the evaluation of the SNR can be found 
in Refs. |llll5j ). As usual, we subject our system to a weak external signal b = bo + A(t) — 
bo + Ab cos(a; s i), rocking the NEP [T5]. In order to have a subthreshold signal, the amplitude 
Ab should satisfy Ab <C b. We have chosen bo as the value of b at which T\<\>^ = T\<$> — 0] = 
when (3 — 0. 

Up to first order in the small amplitude Ab, the transition rates Wi and the functions aj have 
the form indicated in Eqs. (|5j) , ([6[) . But now that depends on the inhomogeneous attractor 
<ps, has the form 

-i 2 



$ = 



1 - 



dx 



After fixing the length of the system L and the kernel range 21 we can use the above indicated 
expressions to find the SNR, that shows the usual bell-shaped form of stochastic resonance as 
a function of the noise intensity. 

In Fig. [9] we show the dependence of the SNR on the kernel range 21 for a fixed value of 2L. 
There is a nonmonotonic behavior in the system's response against variation of 21, that can be 
explained by the following facts: 

— On one hand, the transition rates are decreasing functions of the range 21 for fixed 2L. 
Therefore, the ratio ^1/^2/(^1 + ^2) m the expression for the SNR also reflects this behavior. 
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— The other factor in this expression has a maximum for a kernel range that corresponds to 
"first neighbor" sites, namely x = x' ± I. 



IS. 6 




2 (orb. units) 

Fig. 9. SNR as a function of 21 for 2L fixed. The parameters are D = 0.6, 2L = 6.25, f3 = 0.02 and 
bo = 0.719123. 

The maximum in the system's response as a function of the kernel range is due to the interplay 
between these two factors. From this analysis of the comparative weights of the local (diffusive) 
and nonlocal terms contributing to the SR response, it is apparent that the range of such 
nonlocal kernel has an optimum value yielding a maximum for the SNR 36J . 



5 Conclusions 

We have discussed three different aspects of the phenomenon of stochastic resonance in reaction- 
diffusion systems, within the noncquilibrium potential's framework. In first place we have dis- 
cussed system-size SR in a scalar model. Even though we have not shown the details here, it 
has been also possible to also study other cases In particular, a model of globally coupled 
nonlinear oscillators discussed in [27], showing that it can also be described within the NEP 
framework, with SSSR arising through an "effective" scaling of the noise intensity with the 
system's size. 

In second place we presented a study of SR in systems with a density-dependent (diffusive- 
like) coupling. We initially discuss the case of a scalar system [12], and afterwards extent the 
analysis to an array of FitzHugh-Nagumo units, with a field-dependent activator diffusion |2"2"] . 
For the second system, when both diffusions are constant (that is: D u > and D v — 0), has 
a known form of the NEP [15]. However, in the general case we have not been able to find 
the form of the NEP (but the idea of such a NEP is always underlying our analysis) and have 
to resort to an analysis based on numerical simulations. The result shows that the system's 
response is enhanced due to the particular form of the non-homogeneous coupling. From such 
results, we can conclude that the phenomenon of enhancement of the SNR, due to a selectivity 
in the coupling, initially found for a scalar system |12j is robust, and that the indicated non- 
homogeneous coupling could clearly contribute to enhance the SR phenomenon in very general 
systems. 

Finally, we analyzed an activator-like field including a nonlocal contribution that arise 
through an effective adiabatic elimination of an auxiliary (inhibitor-like) field. By exploiting the 
knowledge of the nonequilibrium potential in such a case, we have analyzed the dependence of 
the SNR on the nonlocal interaction kernel range, founding that there is an optimal value of the 
kernel range, yielding a maximum in the system's response, corresponding to a very localized 
interaction. 

The indicated results clearly show that the "noncquilibrium potential" (even if not known 
in detail [10 ) offers a very useful framework to analyze a wide spectrum of characteristics 
associated to SR in spatially extended or coupled systems. For instance, within this framework, 
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the phenomenon of SSSR looks — as other aspects of SR in extended systems jTT] — as a natural 
consequence of a breaking of the symmetry of the NEP [57] . 



Appendix: Brief review of the nonequilibrium potential scheme 

Loosely speaking, the notion of NEP is an extension to nonequilibrium situations of that of 
equilibrium thermodynamic potential. In order to introduce it, we consider a general system of 
nonlinear stochastic equations (admitting the possibility of multiplicative noises) 

q u = K v (q)+gy(q)^(t), i/=l,...,n; (14) 

where repeated indices are summed over. Equation (Jl4j) is stated in the sense of ltd. The {&(*)}' 
i = 1, . . . , m < n are mutually independent sources of Gaussian white noise with typical strength 
7- 



Graham's approach 



The Fokker-Planck equation corresponding to Eq. ( 14 1 takes the form 

d 2 



7 



2 dq v dqf 



Q^{q)P 



(15) 



where P(q, t; 7) is the probability density of observing q = (qi 7 ...,q n ) at time t for noise 
intensity 7, and Q ufJ "(q) = gi(q)gi(q) is the matrix of transport coefficients of the system, 
which is symmetric and non-negative. In the long time limit (t — > 00), the solution of Eq. (15) 
tends to the stationary distribution P s t(q)- According to [IS], the NEP <P(q) associated to Eq. 
( 15 1 is defined by 

<%) 



- lim 7 lnP st (<7,7). 

7— »o 



(16) 



In other words 



P st (q) d n q=Z{q) exp 



*(g) 

7 



0(7) 



da 



where <P(q) is the NEP of the system and Z(q) is defined as the limit 

1 



In Z(q) = lim 

7^0 



lnP st (g, 7 ) + -<%) 



7 



Here dQ q — d n qj \J G(q) is the invariant volume element in the g-space and G(q) is the deter- 
minant of the contravariant metric tensor (for the Euclidean metric it is G = 1). It was shown 
[T^] that <P(q) is the solution of a Hamilton-Jacobi-like equation (HJE) 

<9<P 1 d@ d$ 

and Z(q) is the solution of a linear first-order partial differential equation depending on <P{q) 
(not shown here). 

Equation ( 16 1 and the normalization condition ensure that <P is bounded from below. Fur- 
thermore, from the separation of the streaming velocity of the probability flow in the steady 
state into conservative and dissipative parts, it follows that 



dt 



dq v 



1 d<P d<P 



i.e. <P is a LF for the dynamics of the system when fluctuations are neglected. Under the 
deterministic dynamics, q v = K"(q), <P decreases monotonically and takes a minimum value on 
attractors. In particular, <P must be constant on all extended attractors (such as limit cycles or 
strange attractors) [16] . 
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Ao's approach 

An alternative way to look into this problem is due to Ao |37j . Let us refer again to the system 
in Eq. ( 14 1 . Following |37j , we introduce now the auxiliary matrix 

H- 1 (q) = S(q) + A(q), 

with S(q) a symmetric matrix while A (a) is an antisymmetric one. H _1 (q) is now used to 
rewrite the initial system as 

H-^q) q = H- X (q) K(q, t) + H^q) £(q, t) = -V*(q, t) + r/(q, f), 

where -V*(q,t) = H" 1 ^ K(q, t), ry(q,t) = H" 1 ^^) and H" 1 ^ = -V«P(q,i) + 
rj(q,t). The new stochastic variables, rj(q,t), fulfill 

(r?(q, i)ry T (q, 0) = 2 S(q) S(t - t') = 2 G- 1 (q)Q(q)[G- 1 (q)] T ( 5(t - t'% 

that imposes a condition on the arbitrary definition of S(q) as we have 

[S(q) + A(q)] Q(q) [S(q) - A(q)] = S(q). 

As shown by Ao, the last equation also implies H(q) + H T (q) = 2Q(q). We also have 
V x [H _1 (q) K(q, <)] = V x [— V#(q, t)\ — 0. From the previous equations, we have in principle 
all the needed conditions to determine H _1 (q), and to obtain from it the potential <£(q, t). The 
interesting feature of this approach is that it resorts neither to P s t(q) nor to the small-noise 
limit, thus being applicable in principle to more general situations. 
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